readr

ERP <- rbind(read.csv("../data/EXP1_pt1.csv"), read.csv("../data/EXP1_pt2.csv")) %>% 
  filter_if(~is.numeric(.), all_vars(!is.infinite(.))) %>%
  filter(time >= -100 & time <= 700) 

ERP %<>% gather(electrode, amplitude, Fp1:T8, factor_key = TRUE) %>%
  mutate(time = as.numeric(as.character(format(round(time, 0), nsmall = 0))),
         Subj = as.factor(Subj),
         Block = as.factor(Block))

names(ERP) <- c("subject", "block", "condition", "time", "electrode", "amplitude")

attaching electrodeLocs

electrodeLocs <- read_delim(loc_path, "\t",
  escape_double = FALSE,
  col_names = c("chanNo","theta","radius", "electrode"),
  trim_ws = TRUE)

ERP %<>% mutate(time = as.factor(time)) %>%
  group_by(subject, block, condition, electrode, time) %>% 
  summarise(amplitude.mean = mean(amplitude),) %>% 
  ungroup() %>% mutate(time = as.numeric(as.character(time)))

ERP_testing <- ERP %>% mutate(time = as.double(time), 
                              amplitude.mean = as.double(amplitude.mean))

electrodeLocs %<>% mutate(radianTheta = pi / 180 * theta,
                          x = radius * sin(radianTheta), y = radius * cos(radianTheta))

ERP %<>% left_join(electrodeLocs, by = "electrode") %>%
  mutate(amplitude.mean = as.numeric(amplitude.mean), 
         electrode = as.factor(electrode)) %>%
  filter(amplitude.mean <= 7 & amplitude.mean >= -7) %>%
  mutate(amplitude.mean = as.numeric(format(round(amplitude.mean, 2), nsmall = 2)))

plotting function: significance

significance <- function(elec = "", blck = "PRS") {

  ERP_test <- filter(ERP_testing, block == blck)
  ERP_testelec <- filter(ERP_test, electrode == elec)
  ERP_testwide <- spread(ERP_test, time, amplitude.mean)

  TERP <- ERP_testwide[(ERP_testwide$electrode == elec), 5:ncol(ERP_testwide)]
  cov.TERP = ERP_testwide[(ERP_testwide$electrode == elec), 1:4]
  cov.TERP = droplevels(cov.TERP)
  design <- model.matrix(~C(subject,sum) + condition, data = cov.TERP)
  design0 <- model.matrix(~C(subject,sum), data = cov.TERP)

  fabh <- erpfatest(TERP, design, design0, nbf = NULL)
  
  return(fabh$significant) }

plotting function: significance_block

significance_block <- function(elec = "", cond = "target") {
  
  ERP_test <- filter(ERP_testing, condition == cond)
  ERP_testelec <- filter(ERP_test, electrode == elec)
  ERP_testwide <- spread(ERP_test, time, amplitude.mean)
  
  TERP <- ERP_testwide[(ERP_testwide$electrode == elec), 5:305]
  cov.TERP = ERP_testwide[(ERP_testwide$electrode == elec), 1:4]
  cov.TERP = droplevels(cov.TERP)
  design <- model.matrix(~C(subject,sum) + block, data = cov.TERP)
  design0 <- model.matrix(~C(subject,sum), data = cov.TERP)

  fabh <- erpfatest(TERP, design, design0, nbf = 6)
  
  return(fabh$significant) }

plotting difference curves with plot.diff

plot.diff <- function(elec = "Cz", labels = FALSE, blck = "PRS") {
    ERP.new <- ERP %>% 
    filter(electrode == elec) %>%
    filter(block == blck)
    
    ERP.match <- ERP.new %>% filter(condition == "control") 
    ERP.mism <- ERP.new %>% filter(condition == "target")
    
    ERP.diff <- ERP.match %>% 
      mutate(amplitude.mean = ERP.mism$amplitude.mean - ERP.match$amplitude.mean,
             condition = "difference")
    
    ERP.new <- bind_rows(ERP.new, ERP.diff)

  my.plot <- ggplot(ERP.new, aes(time, amplitude.mean, color = condition, linetype = condition)) +
      geom_line(stat = "summary") +
      xlim(c(-100, 400)) +
      geom_vline(xintercept = 0, linetype = "dashed" ) +
      geom_hline(yintercept = 0, linetype = "dashed") +
      ggtitle(elec) + 
      scale_linetype_manual(values = c(2, 1, 2), guide = F) +
      scale_color_manual(values = c("steelblue", "black", "indianred3"), 
                         name = "", labels = c("Match", "Difference", "Mismatch")) +
      theme(plot.title = element_text(size = 11, face = "bold")) 

  if (labels == TRUE) { 
    my.plot <- my.plot + 
      labs(x = "Time (ms)", y = expression(paste("Amplitude (", mu, "V)")), colour = "") +
      theme_minimal() } 

  if (labels == FALSE) {
    my.plot <- my.plot +
    labs(x = "", y = "", colour = "") +
    theme_void() }
  
  return(my.plot) }

plotting function bae

bae was actually an acronym, but cannot remember what of to save my life!

bae <- function(elec = "Cz", blck = "PRS", labels = FALSE, significant = TRUE) {
  ERP_gae <- ERP %>% 
    filter(electrode == elec) %>%
    filter(block == blck)
  
  if (blck == "PRS") {
    ERP_CTT = c("steelblue", "indianred3") }
  else {
    ERP_CTT = c("skyblue4", "palegreen3") }

    my.plot <- ggplot(ERP_gae, aes(time, amplitude.mean, color = condition)) +
      geom_line(stat = "summary") +
      xlim(c(-100, 400)) +
      geom_vline(xintercept = 0, linetype = "dashed" ) +
      geom_hline(yintercept = 0, linetype = "dashed") +
      ggtitle(elec) +
      scale_color_manual(values = ERP_CTT, name = "", labels = c("Match", "Mismatch")) +
      theme(plot.title = element_text(size = 11, face = "bold")) 

  if (labels == TRUE) { 
    my.plot <- my.plot + 
      labs(x = "Time (ms)", y = expression(paste("Amplitude (", mu, "V)")), colour = "") +
      theme_minimal() } 

  if (labels == FALSE) {
    my.plot <- my.plot +
    labs(x = "", y = "", colour = "") +
    theme_void() }
  
  if (significant == TRUE) {
    
    spts <- significance(elec, blck) 
    spts1 <- as.data.frame(spts)
    
    if (length(spts) != 0) {
      my.plot = my.plot + geom_point(data = spts1, mapping = aes(x = spts, y = -0.3), size = 3, shape = 20, 
      color = sign_col) } }
  
  return(my.plot) }
bae(elec = "F4", labels = T, significant = F) + transition_reveal(time)

plotting function gae

gae <- function(elec = "Cz", cond = "target", labels = FALSE, significant = TRUE) {
  ERP_gae <- filter(ERP, electrode == elec) %>% filter(time %% 5 == 0) %>%
    filter(condition == cond)
  
  if (cond == "target") {
    ERP_CTT <- c("grey25", "indianred3") }
  else {
    ERP_CTT <- c("steelblue", "skyblue4")
  }

  my.plot <- ggplot(ERP_gae, aes(time, amplitude.mean, colour = block)) +
      geom_line(stat = "summary") +
      xlim(c(-100, 400)) +
      geom_vline(xintercept = 0, linetype = "dashed" ) +
      geom_hline(yintercept = 0, linetype = "dashed") +
      ggtitle(elec) +
      scale_colour_manual(values = ERP_CTT, name = "", 
                          labels = c("Baseline", "Target"))

   if (labels == TRUE) { 
    my.plot <- my.plot + 
      labs(x = "Time (ms)", y = expression(paste("Amplitude (", mu, "V)")), colour = "") +
      theme_minimal() } 

  if (labels == FALSE) {
    my.plot <- my.plot +
    labs(x = "", y = "", colour = "") +
    theme_void() }
  
  if (significant == TRUE) {
    
    spts <- significance_block(elec, cond) 
    spts1 <- as.data.frame(spts)
    
    if (length(spts) != 0) {
      my.plot <- my.plot + 
        geom_point(data = spts1, mapping = aes(x = spts, y = -0.3), size = 3, shape = 20, 
      color = sign_col) } }
  
  return(my.plot) }
## Warning: Removed 780 rows containing non-finite values (stat_summary).

scalp template: theme_topo

theme_topo <- function(base_size = 12)
  {theme_bw(base_size = base_size) %+replace% theme(rect = element_blank(), line = element_blank(), axis.text = element_blank(), axis.title = element_blank())}

circleFun <- function(center = c(0,0), diameter = 1, npoints = 100) {
  r = diameter / 2
  tt <- seq(0,2 * pi, length.out = npoints)
  xx <- center[1] + r * cos(tt)
  yy <- center[2] + r * sin(tt)
  return(data.frame(x = xx, y = yy)) }

headShape <- circleFun(c(0, 0), round(max(electrodeLocs$x)), npoints = 100) 
nose <- data.frame(x = c(-0.075, 0, .075), y=c(.495, .575, .495))

functions int_scalp_plot and int_compare_scalps

int.scalp.plot <- function(TP1, TP2, cond = "control") {
  
  ERP_lme <- ERP %>% 
    filter(condition == cond) %>%
    filter(time >= TP1) %>%
    filter(time <= TP2) %>%
    filter(block == "PRS") %>%
    group_by(subject, block, condition, electrode) %>%
    summarise(Amplitude = mean(amplitude.mean),) %>%
    ungroup %>% filter(condition == cond)
  
  gridRes <- 124
  
  electrodeLocs <- read_delim(loc_path,
                              escape_double = FALSE, 
                              col_names = c("chanNo","theta","radius", "electrode"), 
                              trim_ws = TRUE, 
                              delim = "\t")
  
  electrodeLocs %<>% mutate(radianTheta = pi / 180 * theta,
                            x = radius * sin(radianTheta), 
                            y = radius * cos(radianTheta))

  singleTimepoint <- ERP_lme %>% 
    left_join(electrodeLocs, by = "electrode")
  
  tmpTopo <- with(singleTimepoint, interp(x = x, y = y, z = Amplitude, 
                                          xo = seq(min(x)*2, max(x)*2, 
                                                   length = gridRes), 
                                          yo = seq(min(y)*2, max(y)*2, 
                                                   length = gridRes), 
                                          linear = FALSE, 
                                          extrap = TRUE, 
                                          duplicate = TRUE))

  interpTopo <- data.frame(x = tmpTopo$x, tmpTopo$z)
  names(interpTopo)[1:length(tmpTopo$y)+1] <- tmpTopo$y
  interpTopo <- gather(interpTopo, key = y, value = Amplitude, -x, convert = TRUE)
  interpTopo$incircle <- sqrt(interpTopo$x^2 + interpTopo$y^2) < .7 
  interpTopo <- interpTopo[interpTopo$incircle,] 
  maskRing <- circleFun(diameter = 1.42) 

  ScalpPlotT <- ggplot(interpTopo, aes(x = x, y = y, 
                                       fill = Amplitude)) +
    geom_raster() +
    stat_contour(aes(z = Amplitude, linetype = ..level..<0), 
                 colour = "black", size = 0.8,
                 show.legend = FALSE) +
    theme_topo() +
    scale_fill_viridis(option = "viridis", 
                       limits = c(-0.8, 0.8), 
                       guide = "colorbar", 
                       oob = squish) +
    geom_path(data = maskRing, aes(x, y, z = NULL, fill =NULL), 
              colour = "white", size = 6) +
    geom_point(data = singleTimepoint, aes(x, y), size = 1) +
    geom_path(data = headShape, aes(x, y, z = NULL, fill = NULL), size = 1.5) +
    geom_path(data = nose, aes(x, y, z = NULL, fill = NULL), size = 1.5) +
    coord_fixed() +
    theme(plot.title = element_text(hjust = 0.5, lineheight = 0.5))
  return(ScalpPlotT) }

int.compare.scalps <- function(tp1, tp2) {
  control_map <- int.scalp.plot(TP1 = tp1, TP2 = tp2, cond = "control") 
  target_map <- int.scalp.plot(TP1 = tp1, TP2 = tp2, cond = "target")

final_plot <- annotate_figure(ggarrange(
  target_map + ggtitle("Mismatch"),
  control_map + ggtitle("Match"),
  common.legend = TRUE), 
  top = text_grob(paste("akima amplitude interpolation between" , tp1, " and ", tp2, " ms PSO"), 
                                     face = "bold"))

return(final_plot) }